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Abstract 



We present numerical studies of compressible, decaying turbulence, with and 
without magnetic fields, with initial rms Alfven and Mach numbers ranging up 
to five, and apply the results to the question of the support of star-forming in- 
terstellar clouds of molecular gas. We find that, in ID, magnetized turbulence 
actually decays faster than unmagnetized turbulence. In all the regimes that 
we have studied 3D turbulence — super- Alfvenic, supersonic, sub- Alfvenic, and 
subsonic — the kinetic energy decays as {t — io)"'') with 0.85 < rj < 1.2. We 
compared results from two entirely different algorithms in the unmagnetized 
case, and have performed extensive resolution studies in all cases, reaching 
resolutions of 256^ zones or 350,000 particles. We conclude that the observed 
long lifetimes and supersonic motions in molecular clouds must be due to 
external driving, as undriven turbulence decays far too fast to explain the 
observations. 
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Star-forming clouds of interstellar gas emit strongly in molecular emission lines; tem- 
peratures derived from these lines show that the linewidths greatly exceed the thermal 
sound speed Cg in these clouds. With densities of order n ~ 10^ — 10^ protons per cm^, 
the gas in these clouds can be described by an isothermal equation of state due to the 
efficient radiative cooling allowed by the many low- lying molecular transitions Cloud 
lifetimes are of order 3 x 10'' yr 0, while free-fall gravitational collapse times are only 
iff = (1.4 X 10^ yr)(?7,/10^ cm~'^)~^/^. In the absence of non-thermal support, these clouds 
should collapse and form stars in a small fraction of their observed lifetime. Supersonic hy- 
drodynamical (HD) turbulence is suggested as a support mechanism by the observed broad 
lines, but was dismissed because it would decay in times of order t^. A popular alterna- 
tive has been sub- or trans- Alfvenic magnetohydrodynamical (MHD) turbulence, which was 
thought to decay an order of magnitude more slowly However, analytic estimates and 
computational models suggest that incompressihleMRD turbulence decays as {t—to)~'^, 
with a decay rate 2/3 < rj < 1.0, where to is some characteristic time scale, while incom- 
pressible HD turbulence has been experimentally measured to decay with 1.2 < 77 < 2. 
The difference in decay rates between incompressible HD and MHD turbulence is clearly 
not as large as had been suggested for compressible astrophysical turbulence. In this paper 
we use high-resolution, three-dimensional (3D) simulations to compute the decay rates of 
compressible, homogeneous, isothermal, decaying turbulence with supersonic, sub-Alfenic, 
and super- Alfvenic root-mean-square (rms) initial velocities frms? and show that the decay 
rates in these physical regimes, 0.85 < rj < 1.2, strongly resemble the incompressible results. 

a. Previous Work The decay of the kinetic energy Ek of incompressible HD turbulence 
has been explored in some detail. The energy is predominantly held within low wavenumber 
modes. Dissipation, on the other hand, occurs predominantly from the high wavenumber 
modes. Therefore, the decay rate does not generally depend on the details of the dissipa- 
tion process. Rather, it is controlled by the efficiency of energy transfer from the low to 
high wavenumber modes due to vortex interactions, nonlinear wave interactions, or other 
processes. This leads either to the Kolmogorov decay rate rj = 10/7 for the kinetic energy, 
accompanied by a growth in the external scale of the turbulence C oc t^^'' |10], or to the 



closure model law r] = 1 -|- (s — l)/(s + 3) and C oc t^/Cs+a)^ depending on the injected 
energy spectrum in the low wavenumber region {P{k) oc k'^) and the spatial dimension 



D [|11[]. Hence values for the decay rate r] as low as unity, are predicted, much lower than 
Kolmogorov. However, these results are from studies of spatially free turbulence. When tur- 
bulence is confined 0, a much higher rate is found experimentally: 77 — >■ 2. These confined 
turbulence experiments have been modeled with a mean field theory ||12|| . 

The decay behaviour of incompressible MHD turbulence, in contrast, is controversial. A 
two-dimensional (2D) analysis with constant mean square magnetic potential yields rj = 1, 
a result indeed backed up by numerical studies in 2D [Q. A dimensional analysis in 3D 
assuming magnetic helicity invariance [0] yields rj = 2/3, although low-resolution numerical 
results 1^,0] suggest 77 ~ 1. 



Compressibility introduces an alternative type of complexity |r5|. Nevertheless, the 3D 



decay problem has been modelled numerically |jT6|. Although the evolution of the spectral 
development and spatial structures was explored there, the decay rate was only treated in 
passing, and definitive results are hard to derive. 

Unlike terrestrial turbulence, astrophysical turbulence usually involves full MHD com- 
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pressible flow. Numerical models of one- dimensional (ID), isothermal, compressible, strongly 



magnetized, decaying and forced turbulence have been performed by ^7j, who found decay 
rates rather lower than those discussed above. A recent 3D study |T8| broke new ground 
by examining both weakly and strongly magnetized isothermal, compressible turbulence; 
however the rather low energy decay rates shown there were almost entirely dependent on 



the initial conditions chosen, as explained in ||T9|| . 



b. Numerical Techniques For our HD models of strongly supersonic turbulence we use 
two entirely different HD methods — a second-order, Eulerian, finite-difference code, and a 
smoothed particle hydrodynamics (SPH) code — while for our MHD models we use only the 
finite-difference code. 



This finite-difference code is the well-tested MHD code ZEUS pO[) which uses second- 
order advection, and a consistent transport algorithm for the magnetic fields [Q. It 
resolves shocks using a standard von Neumann artificial viscosity, but otherwise includes 
no explicit viscosity, relying on numerical viscosity to provide dissipation at small scales. 
This should certainly be a reasonable approximation for shock-dominated flows, as most 
dissipation occurs in the shock fronts, where the artificial viscosity dominates in any case. 
The relative simplicity of this Eulerian formulation allows us to perform resolution studies 
showing that our major results are, in fact, independent of the resolution, and thus of the 
strength of numerical viscosity. 



SPH is a particle based approach to solving the HD equations [23|, in which the system is 



represented by an ensemble of particles, each carrying mass, momentum, and fluid properties 
such as pressure, temperature, and internal energy special-purpose processor GRAPE to 
accelerate computation of nearest-neighbor lists [^, allowing us to use as many as 350,000 
particles. 

c. Initial Conditions We chose initial conditions for our models inspired by the popular 
idea that setting up velocity perturbations with an initial power spectrum P{k) oc k" in 
Fourier space similar to that of developed turbulence would be in some way equivalent to 
starting with developed turbulence [|18],|16|. Observing the development of our models, it 
became clear to us that, especially in the supersonic regime, the loss of phase information 
in the power spectrum allows extremely different gas distributions to have the same power 



spectrum. For example, supersonic, HD turbulence has been found in simulations |T6[ to 
have a power spectrum a = —2. However, any single, discontinuous shock wave will also 
have such a power spectrum, as that is simply the Fourier transform of a step function, and 
taking the Fourier transform of many shocks will not change this power law. Nevertheless, 
most distributions with a = —2 do not contain shocks. 

After experimentation, we decided that the quickest way to generate fully developed 
turbulence was with perturbations having a flat power spectrum a = with a cutoff at 
moderate wavenumber (typically fcmax = 8). We set up velocity perturbutions drawn from 
a Gaussian random field fully determined by its power spectrum in Fourier space following 
the standard procedure: for each wavenumber k we randomly select an amplitude from a 
Gaussian distribution centered on zero and with width P{k) = Pq^" with k = |k|, and a 
phase between zero and 27r. We then transform the field back into real space to obtain the 
velocity in each zone. This is done independently for each velocity component. For the SPH 
calculation the velocities defined on the grid are assigned onto individual particles using the 
"cloud-in-cell" scheme P3[. In all of our models we take Cg = 0.1, initial density po = 1, and 



3 



we use a periodic grid with sides L = 2 centered on the origin. These parameter choices 
define our unit system. 

d. One- Dimensional Results To verify our numerical methods, we reproduced the ID, 
MHD results of Gammie & Ostriker Fig. |I]a shows the results of a resolution study 
comparable to their Figure 1, with M = 5, initial uniform field parallel to the x-axis, and 
initial rms Alfven number A = frms/^^A = 1, where v\ = B'^/Airpo. Note that t = 20 in 
our units corresponds to t = 1 in theirs. Aside from a rather faster convergence rate in our 
study, attributable to the details of our choice of initial conditions, we reproduce excellently 
their result: a decrease in wave energy -Ewave = -^k + (-By + BD/^ti by a factor of five in one 
sound-crossing time L/cg. 

We then extended our study by examining the equivalent HD problem, as shown in 
Fig. [l|c, only to find that the decay rate of HD turbulence in ID is significantly slower 
than that of MHD turbulence. This appears to be due to the sweeping up of slower shocks 
by faster ones in the HD case, resulting in pure Bergers turbulence, with linear velocity 
profiles between widely separated shocks exactly as predicted p6[. The result is that there 
are very few dissipative regions, and energy is only lost very slowly. In contrast, multiple 
wave interactions occur in the MHD case producing many dissipative regions and so faster 
dissipation. 

Finally we compared ID models with 256 zone resolution to equivalent 3D models with 
256^ zones. The 3D model loses energy far faster than the ID model in both the HD case 
shown in the thick lines in Fig. and the MHD case shown in Fig. The increased 
number of degrees of freedom available in 3D produces more shocks and interaction regions, 
resulting in increased energy dissipation. 

e. Three- Dimensional Results We next performed resolution studies using ZEUS for 
three different cases with no field, weak field and strong field as described in Fig. 0, and 
summarized in Table |I[ The weak field models have an initial ratio of thermal to magnetic 
pressure /3 = 2, while the strong field models have (3 = 0.08. We ran the same HD model 
with the SPH code to demonstrate that our results are truly independent of the details of the 
viscous dissipation, and so that our lack of an explicit viscosity does not affect our results. 
We also ran two models (R, S) with adiabatic index 7 = 1.4, and an isothermal model (T) 
with initial M = 0.1 to provide a point of direct comparison between our results and those 
for incompressible, Navier-Stokes turbulence. 

The kinetic energy decay curves for the four resolution studies are shown in Fig. ^ For 
each of our runs we performed a least-squares fit to the power-law portion of the kinetic 
energy decay curves, and report the corresponding decay rate rj in Table |. These results 
appear converged at the 5-10% level; it is very reassuring that the different numerical 
methods converge to the same result for the HD case. 

We find that highly compressible, isothermal turbulence (Model D) decays somewhat 
more slowly, with rj = 0.98, than less compressible, adiabatic turbulence (Model R), with 
rj = 1.2, or than incompressible turbulence (Model T), with rj = 1.1 (also see 
Adding magnetic fields (Models L, Q, and S) decreases the decay rate somewhat further in 
the isothermal case to 77 ~ 0.9, with very slight dependence on the field strength or adiabatic 
index. 

/. Conclusions We can draw conclusions for turbulence theory from our models that 
have significant astrophysical implications. What we find remarkable is how closely our 
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results resemble the incompressible results, despite the difference in dissipation mechanisms. 
In incompressible hydrodynamics, kinetic energy is dissipated in vortexes at the smallest 
scales while in supersonic compressible turbulence, kinetic energy is dissipated in shock 
waves, and in MHD turbulence kinetic energy is dissipated in the interactions of nonlinear 
Alfven waves, yet the resulting decay rates differ only slightly. The difference between our 
ID and 3D HD results suggests that somehow the space density of dissipative regions is the 
determining factor, and that in 3D it is somehow quite independent of the detailed physics. 

The clear astrophysical implication of these models is that even strong magnetic fields, 
with the field in equipartition with the kinetic energy, cannot prevent the decay of turbulent 
motions on dynamical timescales far shorter than the observed lifetimes of molecular clouds 



27| . The significant kinetic energy observed in molecular cloud gas must be supplied more 
or less continuously. If turbulence supports molecular clouds against star formation, it must 
be constantly driven, by stellar outfiowsw |2^, photoionization p9[, galactic shear ||30 



or 



some combination of these or other sources. 

Some computations presented here were performed at the Rechenzentrum Garching of 
the MPG. ZEUS was used courtesy of the Laboratory for Computational Astrophysics at 
the NCSA. MDS thanks the DFG for financial support. 
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FIGURES 




FIG. 1. Isothermal, M = 5 models with ZEUS, a) Wave energy decay of MHD models with 
A = 1 and resolutions ranging from 32 to 4096 zones; the 256 zone model is highlighted, b) 
Comparison of the same 256 zone ID MHD model (upper line) to 256^ zone 3D MHD Model Q 
(lower line), c) Kinetic energy decay of HD models with resolutions ranging from 32 (lowest) to 
4096 (highest) zones; the 256 zone model is highlighted, d) Comparison of the same 256 zone ID 
HD model (upper thick line) to 256^ zone HD Model D (lower thick line) , and to 256^ zone MHD 
Model Q (thin line) 
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FIG. 2. 3D resolution studies for M=5, isothermal models. ZEUS models have 32'^ (dotted), 
64^ (short dashed), 128^ (long dashed), or 256'^ (solid) zones, while the SPH models have 7000 
(dotted), 50,000 (short dashed), or 350,000 (solid) particles. Panels show a) HD runs with ZEUS, 
b) HD runs with SPH, c)A = 5 MHD runs with ZEUS, and d) A = l MHD runs with ZEUS. 
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TABLES 



TABLE I. Power law of kinetic energy decay rj (with formal errors from the least squares fits) 
for the 3D models discussed. Initial rms Mach number M and Alfven number A, and adiabatic 
index 7 are given, along with the resolution (in zones per side or thousands of particles) and code 
used. Boldface indicates highest resolution. 
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